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ABSTRACT 

We report on a series of numerical examples and compare several 
algorithms for estimation of coefficients in differential equation models. 
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In this brief note we summarize some of our findings [3] from 
numerical studies on certain aspects of ill-poscdness in inverse or 
parameter estimation problems involving differential equation constraints. 
There is a vast literature (which we shall not attempt to discuss here) on 
a number of questions (e.g., lack of existence and/or uniqueness of 
solutions, lack of continuous dependence of solutions on data) related to 
the estimation of parameters even when the constraining systems are 
algebraic equations or ordinary differential equations. Additional 
difficulties arise when one is attempting to estimate functional (i.e. time 
and/or spatially dependent) coefficients in partial differential equation or 
distributed parameter systems. Here we focus on the role that compactness 
of the admissible parameter or coefficient set plays in such problems. Due 
to the limitations of space, our presentation will be sketchy, with ail but 



the expert reader most likely wishing to consult some of our references for 
further elaboration. 

Consider a general least-squares inverse or parameter estimation 
problem: Minimize J(q) = |Cu(q) -z| z over q€Q AD , C Q, subject to the 

constraints >4(q,u(q)) = F. Here C is an observation operator from the state 
space X to an observation or data space Z, is the admissible subset of 
parameters in the space Q, and the parameter dependent operator A defines 
the dynamics that constrain the problem. In the problems of interest to us, 
A = F represents a partial differential equation (elliptic, parabolic, 
hyperbolic) with parameter functions q which depend on time and/or 
spatial coordinates (or even the state u itself in some nonlinear system 
problems). It is now well-understood (e.g. see [1] for a discussion) that a 
compactness hypothesis (Q^ compact in some Q topology) for plays 
an important theoretical role in both convergence (of approximating 
solutions) and stability (continuity of the estimated parameters with respect 
to the data or observations). Here we wish to demonstrate that this 
compactness also plays an important computational role in such problems. 
To do this, we illustrate the basic ideas with one dimensional elliptic 
systems (so we actually have an ordinary differential equation with 
spatially dependent coefficient). We wish to emphasize, however, that our 
findings are most certainly relevant to problems with more complex system 
dynamics (second order parabolic or hyperbolic equations, or higher order 
equations of elasticity). Indeed we have observed the difficulties and 
phenomena we discuss here in a number of these technically more 
challenging problems. 

Turning to a class of concrete examples, we consider minimization of 

J(q) = fo| u (q) ■ u| 2dx U) 

over qeQ AD CQ = C[0,1] subject to 

D(qDu) = f, u(0) = u( 1 ) = 0. (2) 

Here f is assumed known, D = JL , and u are given observations for 

ox 

u = u(q) the solution of (2). For computational purposes, we replace this 
original problem by a sequence of approximating problems where the states 
u are replaced by Galerkin approximations u N (for example, here we take 



u N in the N-l dimensional space of linear splines with grid size 1/N and 
satisfying the boundary conditions u N (0) = u N (l) = 0) and the approximate 
parameters q M are chosen from an approximating set Q M for Q^. That is, 
our algorithms are used to seek q M € Q M C Q = C[0,1] that minimizes 

j N (q M ) = Jo | u N (q M ) - u| 2 dx. (3) 

A convergence theory (as the dimensions of the approximating spline 
spaces increase i.e., N -» « M -♦ °») can be given where one may use either 
linear or cubic splines for the state approximations and for the parameter 
approximations (e.g. see [2], [4] for the ideas). An essential feature of these 
particular convergence proofs is that the admissible parameter set and 

its approximations Q M lie in some compact subset of C[0,1]. This same 
compactness assumption plays a fundamental role in proving stability (e.g., 
continuity of the inverse of the mapping from the parameter estimates to 
the observations or data) as is discussed in [1], for example. 

Perhaps the most direct way to interpret the compactness 

requirements is in terms of constraints on the parameters. For example, in 

the computations reported on herein, we imposed compactness in C[0,1] by 
putting pointwise upper and lower bounds on the parameter function 
values as well as an upper bound on the absolute values of the slope of 
the functions. In practice it is common to ignore functional constraints, 

imposing the pointwise upper and lower bounds to insure that the 

optimization algorithms perform satisfactorily. The results summarized in 
this note illustrate the apparent necessity in many examples of including 
the full compactness constraint in computational algorithms. Examples are 
given here and in [3] where both stability and convergence properties are 
as expected whenever a constrained estimation procedure is employed 
whereas instability and divergence are in evidence when unconstrained 
techniques are used. It is safe to speculate that similar behavior occurs in 
problems with parabolic and hyperbolic as well as elliptic systems. In our 
own work and in that reported in the literature - e.g. see Yoon and Yeh 
[6], one sometimes encounters severe problems with oscillations in the 
estimates for q as one pushes the algorithms for increased accuracy in the 
parameter estimates (i.e. as one lets M -* °°). As the examples in this note 
and [3] demonstrate, these difficulties can to* some extent be alleviated by 
imposition of compactness constraints. 



An alternative but essentially theoretically equivalent approach 
involves the use of Tikhonov regularization as formulated by Kravaris and 
Seinfeld in [5]. One restricts the parameter set to Q R C Q with Q R 
compactly imbedded in Q and then modifies the original least squares 
criterion J to minimize Jg = J + 0jq| R where |-| R is the norm in Q R and 
0 is a regularization parameter. Thus minimizing sequences for Jg are 
bounded in Q R and hence compact in Q; this is, in some sense, roughly 
equivalent to minimizing J over a restriction of Q which is compact even 
though the minimization of Jg only produces (hopefully) an approximation 
to the minimizer for the original criterion J. In the cases considered below, 
we use Q r = H 1 while Q = C (which corresponds to A = C 1 and R = H 2 
in the notation of [5]). 

As we shall see below, each approach has inherent difficulties in 
choosing related imbedding parameters: in the first, the estimates produced 
are sensitive to the constraints (the bound L on the derivatives of the 
parameters in the computations summarized here) while the estimates 
produced using regularization are quite sensitive to the regularization 
parameter 0. 

We carried out a series of numerical tests to compare spline based 
algorithms (linear spline approximations for both the states and parameters) 
on a number of examples for three cases: the unconstrained minimization 
of J; the constrained minimization of J; and unconstrained minimization of 
a regularized criterion Jg. Details of our packages and the algorithms are 
given in [3]. Here we only note that for the constrained minimization we 
used a reduced gradient algorithm with a corresponding gradient algorithm 
for the unconstrained minimization. For the compactness constraints we 
used |Dq M (x)| ? L and .5 5 q M (x) S 10.0 in all our examples. 

Our algorithms were compared on examples for which we knew the 

* A 

true solutions, i.e. we used "true" parameter values q to generate data u 
(in some cases with noise) as described in [3]. We summarize our findings 
and present representative results. 

CONVERGENCE 

In Figure 1 we compare estimates for several values of approximation 
indices N, M produced for an example (Example 2 of [3]) with 
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CONSTRAINED (L=l) UNCONSTRAINED 



(N,M) = (8,15) 




(N,M) = (16,16) 



(N,M) = (64,32) 


FIGURE 1 








FIGURE 2 




TIKHONOV ESTIMATES AS 0 VARIES 


CONSTRAINED ESTIMATES AS L VARIES 
(N,M) = (64,16) 



FIGURE 3 
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q*M 


2 + x 
8/3 - x 
4/3 + x 


0 S x « 1/3 
1/3 S x S 2/3 
2/3 <= x S 1, 


(4) 


u(x) = i/x(l - vOC) , 


(5) 


in the system (2). We chose L = 1 (the exact maximum of |Dq*(x)|) for 
the constrained algorithms. For low values of N (often a desirable situation 
in practice) compared to M, the unconstrained estimates are totally useless. 
For larger values of N, the unconstrained estimates are improved with only 
small oscillations appearing at each end. In all cases, convergence took 
much longer for the unconstrained package. For this same example we 
depict in Figure 2 the estimates obtained using Tikhonov regularization 
with N = 64, M = 16 and several values of the regularization parameter 6. 
Results for a slightly different example (Example 3 of [3]) with the same 
u but q*(x) piecewise linear as in (4) except with slopes T 2 are depicted 
in Figure 3. Here we illustrate, for N = 64, M = 16, the typical 

performance of Tikhonov regularization as 0 changes and that of the 
constrained minimization as L varies. As one might expect, the estimates 
begin to resemble unconstrained estimates as 3 -* 0 and L -» °». 

STABILITY 

To investigate stability with respect to noise in the data, we took the 
* 

true u, q and f associated with (4), (5) above but perturbed u to produce 
data u = u p for the least squares criterion J N of (3). We used perturbations 

of the form u p (x) = u(x) + p(x)/K where p(x) is a perturbing function and 

K can be varied to control the size of the perturbation. As R ■* ”, we 
have u p (x) -» u(x) for bounded perturbations p. In our numerical 
experiments we used two different perturbation functions: p t (x) = 

x(l - x), p 2 (x) = 1. Note that p t (like u) satisfies the homogeneous 

boundary conditions while p 2 does not. The unconstrained, constrained and 
^Tikhonov estimates for several values of K with perturbation function Pj 
in the data and N = 8, M = 15 are given in Figure 4. The depicted 

behavior is typical: For all values of N the behavior of the constrained 

and Tikhonov methods are similar, with the estimates improving steadily as 
u p (x) -* u(x), i.e. as the noise in the observations tends to zero. In Figure 
5, we present results obtained using the perturbation p 2 with the 
unconstrained, constrained and Tikhonov estimation procedures for several 
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UNCONSTRAINED ESTIMATES 


K=IOO 



TIKHONOV ESTIMATES 


K=100 



CONSTRAINED ESTIMATES 


(8,15) , PERTURBATION p 
FIGURE 4 









values of (N,M) : (8,15), (16,16) and (64,16). The L” norm of the error 
(from the true values q*) in the final estimates is graphed versus the 
values of K. 

SUMMARY REMARKS 

The results here and in [3] demonstrate severe problems in some 
instances with using an unconstrained algorithm to estimate the parameter 
q in examples such as (1), (2). When modified, either by regularizing the 
problem using Tikhonov regularization or by constraining the estimate set 
as in this note, the algorithm does give good estimates. 

Unlike the unconstrained algorithm, both the Tikhonov and 
constrained algorithms are stable with respect to increasing M while 
holding N fixed. However as N is increased the estimates from the 
Tikhonov algorithm do not improve as much as do those of the constrained 
algorithm. The Tikhonov estimates are biased by the regularization of the 
cost functional, and never show all the detail of q when q has significant 
variation. 

Both the constrained and Tikhonov estimation algorithms are stable 
with respect to systematic errors in the observation data, while, except 
when N is large, the unconstrained algorithm fails to give good results on 
even the exact data. 

For both the Tikhonov and constrained algorithms there are 
parameters which affect the algorithm’s performance. For the constrained 
algorithm suitable constraints must be found while for the Tikhonov 
algorithm suitable values of 6 must be found. The constrained algorithm 
has the advantage that the constraints used here, i.e. limits on the slope of 
q, have an obvious meaning, and so may well be (at least approximately) 
known in advance. In the Tikhonov algorithm 0 has no obvious meaning. 
It must be chosen by looking at the change in the estimate behavior as 0 
changes and perhaps using some a priori knowledge about the shape of q 
to choose values of 0 that give an estimate that is neither too flat, nor 
too oscillatory. For the constrained algorithms, the estimates are sensitive to 
the slope constraint parameter L. We have begun investigations into how 
one might use this sensitivity in some type of adaptive manner in 
algorithms to choose a "best" value of L (and hence a good parameter 
estimate). In Figure 6 we depict some of our initial findings. In this figure 










we graph the square of the H 1 norm of the final estimate versus the 

constraint L for several different examples and various values of N (with 

M fixed at 16). All but the second of the examples involve true 

parameters q* of the form (4), differing only in the slope of the piecewise 

linear functions. In the first example Dq*(x) = T 1, while in the last two 
* * 

Dq (x) = T 2 and Dq (x) = T 5 respectively. The second is made up of 

£ 

piecewise linear and parabolic segments satisfying |Dq (x)| $ 1. Note that 
it is not necessary to know the true values q in order to obtain the 

graphs in this figure. Furthermore, we observe a striking separation in the 
values of the H 1 norm of the estimate for q M at the value of L 
corresponding to the desired value of L to be used with each example. We 

are continuing our investigations into how these and other features of 

some of our results might be used to develop "adaptive" constrained 
parameter estimation algorithms. 
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